
function [dfdt] = waves(t,y)
global Nx Ny L C E F;
dfdt(1:2*Nx*Ny,1) = -L*y(1:2*Nx*Ny,1) ;%+ C * y(2*Nx*Ny+1:3*Nx*Ny,1);
%F = C(1:Nx*Ny,1:Nx*Ny);
%dfdt(2*Nx*Ny+1:3*Nx*Ny,1) = E*y(2*Nx*Ny+1:3*Nx*Ny,1) + F*y(Nx*Ny+1:2*Nx*Ny,1);